/*
 ============================================================================
 Name        : MC-Pi.c
 Author      : Aleksandar Stancic
 Version     :
 Copyright   :
 Description : Hello World in C, Ansi-style
 ============================================================================
 */

//#define MPI_ENABLED

#define DEFAULT_NUM (1000)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#ifdef MPI_ENABLED
	#include <mpi/mpi.h>
#endif
#include <limits.h>

unsigned long realEvals = 0;
unsigned long buckets[100] = { 0 };

#define PI 3.1415926535897932384626433832795028841971693993751058209749

long double evalFunction(long double x)
{
	realEvals++;
	//buckets[(int)(x * 100)]++;
	//return sqrt(1.0 - x * x);
	//return sin(1.0 / x);
	//printf("Samp X = %2.5Lf ; sin(x) = %2.5lf\n", x, sin(x));
	return sin(x);
	//return x * x * cos(x);
	//return 1.0 / (1.0 + x);
}

double func3d (long double x, long double y, long double z)
{
    long double A = 1.0 / (PI * PI * PI);
    return A / (1.0 - cos (x) * cos (y) * cos (z));
}

void integrate(long double x1, long double x2, unsigned long samples,
		long double *avg, long double *variance)
{
	unsigned long i;
	long double sum = 0.0;
	long double multi = (x2 - x1) / (long double)RAND_MAX;
	long double diff = (x2 - x1);
	long double last = evalFunction(rand() * multi + x1);
	long double current;
	long double var = 0.0;
	long double totVar = 0.0;
	long double sx = rand() * multi;
	int swch = 0;
	sum = last;
	printf("Max samples : %ld\n", samples);
	for(i = 0; i < samples - 1; i++)
	{
		//if(swch)
		//	sx = diff - sx;
		//else
			sx = rand() * multi;
		//swch = !swch;
			if(samples > (unsigned long)100) return;
		current = evalFunction(sx + x1);
		sum += current;
		var = current - last;
		totVar += var * var;
		printf("Samp : %ld\n", i);
		//var += current * current;
	}
	printf("totVar : %Lf\n", totVar);
	totVar /= samples;
	sum /= samples;
	sum *= diff;
	//totVar = var / (samples - 1) - sum * sum;


	*avg = sum;
	*variance = totVar;
}

long double uniformSampling(long double x1, long double x2, unsigned long samples)
{
	long double sum = 0.0;
	for(unsigned long i = 0; i < samples; i++)
		sum += evalFunction(x1 + (x2 - x1) * ((long double)i / (long double)(samples - 1)));
	printf("Uniform sum : %3.5Lf\n", sum);
	sum /= (long double) samples;
	sum *= (x2 - x1);
	return sum;
}

long double random3d(long double x1, long double y1, long double z1,
		long double x2, long double y2, long double z2, unsigned long samples)
{
	long double sum = 0.0;
	long double multiX = (x2 - x1) / (long double)RAND_MAX;
	long double multiY = (y2 - y1) / (long double)RAND_MAX;
	long double multiZ = (z2 - z1) / (long double)RAND_MAX;
	//long double volume = (x2 - x1) * (y2 - y1) * (z2 - z1) *
	//for(unsigned long i = 0; i < samples; i++)
		//sum += func3d(rand() * multiX + x1, rand() * multiY + y1, rand() * multiZ + z1) *
		//;
	return sum / (long double)samples;
}

long double recIntegrate(long double x1, long double x2, unsigned long samples,
		long double fraction, int step)
{
	if(step == 0 || samples < 2)
	{
		long double avg, var;
		//printf("Step 0\n");
		integrate(x1, x2, samples, &avg, &var);
		//printf("Step 0 over\n");
		return avg * fraction;
		//return avg;
	}
	else
	{
		long double avg1, var1, avg2, var2, avg3, avg4, vsum, nres;
		//printf("Int range : %.5lf - %.5lf\n", x1, x2);
		fraction *= 0.5;
		integrate(x1, (x1 + x2) / 2.0, (long)ceil(samples / 4.0), &avg1, &var1);
		integrate((x1 + x2) / 2.0, x2, (long)ceil(samples / 4.0), &avg2, &var2);
		vsum = var1 + var2;

		/*
		if(var1 > var2)
		{
			nres = recIntegrate(x1, (x1 + x2) / 2.0, samples / 2, fraction * 0.5, step - 1);
			return ((avg1 + nres) / 2.0 + avg2) / 2.0;
		}
		else
		{
			nres = recIntegrate((x1 + x2) / 2.0, x2, samples / 2, fraction * 0.5, step - 1);
			return (avg1 + (avg2 + nres) / 2.0) / 2.0;
		}
		 */


		if(vsum == 0.0)
		{
			var1 = 0.5;
			var2 = 0.5;
		}
		else
		{
			var1 = var1 / vsum;
			var2 = var2 / vsum;
		}
		printf("Var1 : %5.3lf    Var2 : %5.3lf    Fraction : %5.5lf\n", var1, var2, fraction);
		nres = recIntegrate(x1, (x1 + x2) / 2.0, (long)ceil(samples / 4.0 * var1), fraction, step - 1); // * var1;
		nres += recIntegrate((x1 + x2) / 2.0, x2, (long)ceil(samples / 4.0 * var2), fraction, step - 1); // * var2;
		//return ((avg1 + avg2) / 2.0 + nres) / 2.0; //((avg1 + avg3) / 2.0 + (avg2 + avg4) / 2.0) / 2.0;
		//printf("First : %3.5Lf    Second : %3.5Lf\n", (avg1 + avg2) * fraction, nres);
		return ((avg1 + avg2) * fraction + nres) / 2.0;

	}
}

int main(int argc, char* argv[]) {
	unsigned long max = DEFAULT_NUM;
	long double accum = 0, total = 0, var;
	int rank = 0, size = 1;
	int steps = 4;
	int uniform = 0;
	int dim3 = 0;
	int random = 0;

	int timeStart = 0;
	int fmode = 0;
	int i;

	for(i = 1; i < argc; i++)
	{
		if(strcmp(argv[i], "-max") == 0 && i + 1 < argc)
		{
			max = atol(argv[i + 1]);
			if(max < 1)
			{
				printf("Broj iteracija mora biti veci od nule. (Default je %d)\n", DEFAULT_NUM);
				return 1;
			}
		}
		if(strcmp(argv[i], "-steps") == 0 && i + 1 < argc)
		{
			steps = atoi(argv[i + 1]);
			if(steps < 1)
			{
				printf("Broj koraka mora biti veci od nule. (Default je %d)\n", DEFAULT_NUM);
				return 1;
			}
		}
		else if(strcmp(argv[i], "-lmax") == 0)
		{
			max = ULONG_MAX;
			printf("ULONG_MAX na lokalnoj masini : %lu\n", ULONG_MAX);
		}
		else if(strcmp(argv[i], "-float") == 0)
		{
			fmode = 1;
		}
		else if(strcmp(argv[i], "-uniform") == 0)
		{
			uniform = 1;
		}
		else if(strcmp(argv[i], "-random") == 0)
		{
			random = 1;
		}
		else if(strcmp(argv[i], "-3d") == 0)
		{
			dim3 = 1;
		}
	}

#ifdef MPI_ENABLED
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	MPI_Comm_size(MPI_COMM_WORLD, &size);
#else
	timeStart = time(0);
#endif

	srand(time(0) * (rank + 1));

	if(dim3)
	{
		accum = random3d(0, 0, 0, PI, PI, PI, max / size) * 8.0 * PI;
	}
	else if(random)
	{
		//accum = random3d(0, 0, 0, PI, PI, PI, max / size) * 8.0;
		integrate(0.0, PI, max / size, &accum, &var);
	}
	else if(uniform)
	{
		accum = uniformSampling(0.0, PI, max / size);
	}
	else
	{
		printf("Start\n");
		accum = recIntegrate(0.0, PI, max / size, PI, steps);
		printf("End\n");
	}


#ifdef MPI_ENABLED
	MPI_Reduce(&accum, &total, 1, MPI_LONG_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
	accum /= size;
#else
	//count = localCount;
#endif

	if(rank == 0)
	{
		max *= size;
		printf("%.15Lf\n", accum);
		printf("Samples : %ld\n", realEvals);

		//for(int i = 0; i < 100; i++)
		//	printf("%3d - %3d : %10ld\n", i, i + 1, buckets[i]);
	}

#ifdef MPI_ENABLED
	MPI_Finalize();
#else
	printf("Time elapsed : %ld s\n", time(0) - timeStart);
#endif

	return EXIT_SUCCESS;
}
